Multi-omics data integration

# options
options(stringsAsFactors = F)
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
knitr::opts_knit$set(progress = FALSE)

We feed log transformed values, using log1p() base R function, for accessibility at 15,000 regions in the genome with the most variance in chromatin accessibility, all transcripts and all proteins with abundance measurements into MOFA+ for model generation.

Figure 5A: MOFA factors overview

# create the MOFA object
all_df_shared_MOFAobject <- create_mofa(all_df_shared)

# Load in model results generated using _src/MOFA_train_model_CPU.r
all_df_shared_model <- load_model(here("../pQTL_website/_data/all_df_shared_30factors_2022-04-11.hdf5"), remove_inactive_factors = F)
# Warning message:
# In .quality_control(object, verbose = verbose) :
#   Factor(s) 3 are strongly correlated with the total number of expressed features for at least one of your omics. Such factors appear when there are differences in the total 'levels' between your samples, *sometimes* because of poor normalisation in the preprocessing steps.

# Removing Factor 3 that shows high correlation to the total # of expressed features.
all_df_shared_model <- subset_factors( all_df_shared_model, factors = c(1:2, 4:30))

# Removing Factors that don't explain at least 1% of variation in a data set. 
all_df_shared_factor_var <- all_df_shared_model@cache$variance_explained$r2_per_factor %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(single_group.Chromatin > 1 |
           single_group.Protein > 1 |
           single_group.Transcript > 1) %>%
  column_to_rownames()
all_df_shared_model_filtered <- subset_factors(all_df_shared_model,
                                        factors = as.numeric(gsub("Factor","",rownames(all_df_shared_factor_var))))

# updating the sexes in the metadata to use Female/Male instead of 0s and 1s.
merged_metadata <- merged.covar2 %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  select(-lifr) %>%
  rename("sample"="rowname")
shared_samples_metadata <- all_df_shared_model_filtered@samples_metadata %>%
  left_join(.,merged_metadata) %>%
  mutate(sex = ifelse(sex ==0,"Female", "Male"))
samples_metadata(all_df_shared_model_filtered) <- shared_samples_metadata


# Get % variance explained for each factor
all_df_shared_var_explained <- (calculate_variance_explained(all_df_shared_model_filtered))

# Get Factor weights
all_df_shared_factors <- get_factors(all_df_shared_model_filtered,
  factors = "all",
  as.data.frame = T
)
# Add sample details and convert to matrix with Factors in rows and samples in columns.
all_df_shared_factors_mat <- all_df_shared_factors %>%
  pivot_wider(id_cols = "sample", names_from = "factor", values_from = "value") %>%
  left_join(., select(threeway.shared.samples, top_muga, sampleid), by =c("sample"="top_muga")) %>%
  column_to_rownames("sample") %>%
  select(-sampleid) %>%
  as.matrix()

# Get feature weights
all_df_shared_weights <- get_weights(all_df_shared_model_filtered, as.data.frame = TRUE, scale = TRUE)

# Correlate Factors with covariates and get R values
correlate_factors_with_covariates(all_df_shared_model_filtered,
  covariates = c("sex","lifr_geno"),
  plot="r",
  return_data = TRUE) %>%
  as_tibble( rownames = "Factor") %>%
  pivot_longer( cols = c("sex","lifr_geno"),
                names_to = "covariate",
                values_to = "pearson_r") -> corr_to_cov_r

# Correlate Factors with covariates and get p-values
correlate_factors_with_covariates(all_df_shared_model_filtered,
  covariates = c("sex","lifr_geno"),
  plot="log_pval",
  return_data = TRUE) %>%
  as_tibble( rownames = "Factor") %>%
  pivot_longer( cols = c("sex","lifr_geno"),
                names_to = "covariate",
                values_to = "log_pval") -> corr_to_cov_logpval

# merge R & p-values into a single data frame
corr_to_cov_r %>%
  full_join( corr_to_cov_logpval) -> all_df_shared_corr_covar


# get atac-seq peak drivers + run overrepresentation analysis using LOLA

# these are all the ATAC-seq peaks that are fed into MOFA that we will use as the custom background in overrepresentation analysis
background_atac_peaks <-  all_df_shared %>%
  filter(view == c("Chromatin"))  %>%
  rename( peak_id = feature) %>%
  select( -view,-value, -sample) %>%
  distinct() %>%
  separate( peak_id, into = c("Chr", "Start","End"), remove = FALSE) %>%
  mutate( Chr = gsub("peak","chr",Chr)) %>%
  makeGRangesFromDataFrame(.,
                               keep.extra.columns = F,
                               seqnames.field = c("Chr"),
                               start.field = "Start",
                               end.field = "End")

get_tf_ora <- function( factor_df, factor_num, bg_peaks){
  
  # get top ATAC-seq peak drivers
  top_atac_drivers <- factor_df %>%
    slice( which( factor_df$value %in%
        boxplot.stats( filter(factor_df, factor == factor_num & view =="Chromatin")$value)$out )
    )
  
  # convert the top ATAC-seq peak drivers to a compatible object for LOLA
  atac_peaks_for_lola <- top_atac_drivers %>%
    rename( peak_id = feature) %>%
    separate( peak_id, into = c("Chr", "Start","End"), remove = FALSE) %>%
    mutate( Chr = gsub("peak","chr",Chr)) %>%
    makeGRangesFromDataFrame(.,
                               keep.extra.columns = F,
                               seqnames.field = c("Chr"),
                               start.field = "Start",
                               end.field = "End")
  # run LOLA with custom background
  lola_results <- runLOLA(atac_peaks_for_lola,
                          bg_peaks,
                          regionDB,
                          cores=1)
  # add q-value for filtering later
  lola_results$qValue <- (qvalue( 10^(-lola_results$pValueLog )))$qvalues

  return(lola_results)
  
}

# get_tf_ora( all_df_shared_weights, "Factor1", background_atac_peaks)
# I ran the function above for all factors and saved the output
# below is the filtered results for qvalue <0.05 and cell type as ESCs
load( here("../pQTL_website/_data/","MOFA_Factor_LOLA_results.RData")) # lola_combined

Figure5A_data_var_exp <- all_df_shared_var_explained$r2_per_factor$single_group %>%
  as_tibble(rownames = "Factor") %>%
  mutate(factor_num = as.numeric(gsub( "Factor" ,"",Factor)) )%>%
  arrange((factor_num)) %>%
  pivot_longer(c("Chromatin","Protein","Transcript")) %>%
  mutate( Factor = factor(Factor, levels = unique(Factor))) %>%
  mutate( name = factor( name, levels = c("Chromatin","Transcript","Protein"))) %>%
  select( `Data set` = name, Factor, `% Variation Explained` = value)

Figure5A_data_cor_to_covar <- all_df_shared_corr_covar %>%
  mutate( Factor = factor(Factor, levels = unique(Factor))) %>%
  mutate( covariate = ifelse( covariate =="sex", "Sex", "LIFR genotype")) %>%
  mutate( covariate = factor( covariate, levels = c("Sex", "LIFR genotype"))) %>%
  select(Factor, covariate, Correlation =pearson_r)

Figure5A_data_tf_ora <- lola_combined %>% 
  filter( antibody %in% c("Nanog","Sox2","Pou5f1")) %>% 
  filter( description %in% c("ES cells expressing control shRNA targeting GFP",
                             "V6.5 (C57BL/6-129) murine ES cells were grown under typical ES conditions on irr",
                             "N/A")) %>% 
  filter( !description %in% c("ES cells expressing shRNA targeting Kdm4b",
                             "ES cells expressing shRNA targeting Kdm4c",
                             "ES cells cultured in FGF4; heparin and Activin A for 1 day after Oct3/4 deleted",
                             "ES cells differentiated in vitro for 2 days; Brn2 expression induced using Tet-O",
                             "ES cells cultured in FGF4; heparin and Activin A for 2 days after Oct3/4 deleted",
                             "ES cells cultured in FGF4; heparin and Activin A for 3 days after Oct3/4 deleted",
                             "ES cells differentiated in vitro for 2 days")) %>% 
  full_join( tibble(
                factor = c( paste0("Factor ", seq(1:23))),
                factor_num = seq(1:23)) 
              ) %>% 
  group_by(factor, antibody, factor_num) %>% 
  dplyr::summarize( mean_odds = mean(oddsRatio, na.rm= T)) %>% 
  ungroup() %>% 
  arrange(factor_num) %>% 
  mutate( Factor = factor(factor, levels = unique(factor))) %>%
  select(Factor, TF = antibody, `Odds ratio` = mean_odds)

# % variation explained heatmap.
Figure5A_data_var_exp %>% 
  ggplot()+
  aes( 
       y = Factor,
       x = `Data set`,
       fill = `% Variation Explained`)+
  geom_tile()+
  xlab("")+
  ylab("")+
  theme_pubclean( base_size = 18)+
  scale_fill_gradient2( limits = c(0, 10))+
  #scale_size_area()+
  labs(fill = "% Variation Explained", col = "")+
  scale_x_discrete( expand = expansion(mult = 0))+
  theme(
        legend.position = "top", 
        legend.title = element_text(size = 18, vjust = 0.9),
        legend.text  = element_text(angle=45, size =12, vjust = 0.5),
        )-> all_df_shared_var_plot

# correlation to covariates.
Figure5A_data_cor_to_covar %>%
  filter(covariate == "Sex") %>% 
  ggplot()+
  aes(
    x = covariate,
    y = Factor,
    fill = abs(Correlation)
  )+
  geom_tile()+
  xlab("")+
  ylab("")+
  theme_pubclean(base_size = 18)+
  scale_fill_gradient(low = "white",high = "dark red", limits = c(0,1))+
  theme(legend.position = "top",
        axis.ticks.y = element_blank(),
        legend.text  = element_text(angle=45, size =12, vjust = 0.5),
        legend.title = element_text(size = 18, vjust = 0.9),
        )+
  scale_y_discrete( labels = NULL)+
  scale_x_discrete( expand = expansion(mult = 0))+
  labs( fill = "Correlation") -> all_df_shared_covar_plot_sex

Figure5A_data_cor_to_covar %>%
  filter(covariate != "Sex") %>% 
  ggplot()+
  aes(
    x = covariate,
    y = Factor,
    fill = abs(Correlation)
  )+
  geom_tile()+
  xlab("")+
  ylab("")+
  theme_pubclean(base_size = 18)+
  scale_fill_gradient(low = "white",high = "deeppink4", limits = c(0,1))+
  theme(legend.position = "top",
        axis.ticks.y = element_blank(),
        legend.text  = element_text(angle=45, size =12, vjust = 0.5),
        legend.title = element_text(size =12)
        )+
  scale_y_discrete( labels = NULL)+
  scale_x_discrete( expand = expansion(mult = 0))+
  labs( fill = "") -> all_df_shared_covar_plot_lifr

# Odds ratio for TF binding overrepresentation for NANOG, OCT4, SOX2 in ATACseq peak drivers
Figure5A_data_tf_ora %>% 
  mutate( mean_odds_capped = ifelse( `Odds ratio` >10, 10, `Odds ratio`)) %>% 
  ggplot()+
  aes(
    x = TF,
    y = Factor,
    fill = mean_odds_capped
  )+
  geom_tile()+
  xlab("")+
  ylab("")+
  theme_pubclean(base_size = 18)+
  scale_fill_gradient(low = "light gray",high = "#ffa08aff", limits = c(0,10)) +
  theme(legend.position = "top",
        legend.text  = element_text(angle=45, size =12, vjust = 0.5),
        legend.title = element_text(size = 18, vjust = 0.9),
        axis.text.x =  element_text(angle = 0),
         axis.ticks.y = element_blank(),
        )+
  scale_y_discrete( labels = NULL)+
  scale_x_discrete( expand = expansion(mult = 0))+
  labs( fill = "Odds ratio") -> lola_sum_figure

figure5a <- ggarrange( all_df_shared_var_plot,
                       all_df_shared_covar_plot_sex,
                       all_df_shared_covar_plot_lifr,
                       lola_sum_figure,
                       widths = c(0.45,0.15,0.15,0.2), nrow = 1,
                       labels = "A",
                       font.label = list( size = 28))


figure5a

MOFA yielded 23 latent factors that capture variation in one or more layers of genomic data. For each factor, percent of variation explained in chromatin accessibility, transcript abundance, and protein abundance is displayed as a heatmap, as is the correlation of each factor to experimental covariates including sex and genotype at the Lifr locus. Heatmap on the right indicates overrepresentation of pluripotency regulator binding sites (NANOG, OCT4 (Pou5f1) and SOX2) among the top chromatin drivers of each factor.

Figure5A_data_var_exp %>% 
  mutate_if( is.numeric, round, 2) %>% 
  create_dt()

Figure5A_data_cor_to_covar %>% 
  mutate_if( is.numeric, round, 2) %>% 
  create_dt()

Figure5A_data_tf_ora %>% 
  mutate_if( is.numeric, round, 2) %>% 
  create_dt()

Table S7: MOFA factor values for samples and lists of molecular features with their weights for each of 23 MOFA factors.

MOFA factor values (columns) for each of the 163 samples (rows) are listed in the first tab. Weights for each of the 23 MOFA factors (columns) are listed for every protein, transcript, and chromatin peak (rows), with tabs corresponding to the three data types.Feature annotations include Ensembl protein identifier and gene symbol for proteins, Ensembl gene identifier and gene symbol for transcripts, and chromatin peak identifier and proximal gene annotation from R/ChIPseeker for ATAC-seq peaks.

Download Table_S7.xlsx

Figure S5A: Drivers of MOFA factors.

all_df_shared_weights %>% 
  filter( abs(value) > 0.01) %>%
  #group_by(view) %>% summarise( value = median(value))
  group_by(factor) %>% 
  count( view, .drop = FALSE) %>%
  ggplot()+
  aes(
    x = view,
    y = factor,
    fill = n
    )+
  geom_tile()+
  xlab("")+
  ylab("")+
  theme_pubclean( base_size = 18)+
  #scale_color_viridis_c(option = "C")+
  scale_fill_gradient2(limits = c(0, 15000))+
  labs(fill = "# of features", col = "")+
  scale_x_discrete( expand = expansion(mult = 0))+
  theme(
        legend.position = "top", #c(0.75,0.95),
        #axis.line.y = element_blank(),
        legend.title = element_text(size = 18, vjust = 0.9),
        legend.text  = element_text(angle=45, size =12, vjust = 0.5),
        ) 

Figure S5. Heatmaps showing the number of features in each data set with abs(weight) > 0.01 for 23 MOFA Factors.

Figure 5B: QTL mapping with MOFA Factors

# load from file
knitr::include_graphics(here("Figure5B.png"))

Figure 5B: Above shows a depiction of QTL mapping with MOFA factors to identify the genetic modifiers driving variation across three molecular layers. Below is a table of QTL peaks that map above the genome-wide significance threshold calculated individually for each factor. Loci that were previously observed as molecular QTL hotspots are denoted in the “Type” column.

# rankZ transform before mapping
all_df_shared_factors_mat_rankZ <- apply( all_df_shared_factors_mat, 2, rankZ)



# run qtl scans with sex as an additive covariate
all_df_scans <- scan1( genoprobs = shared.probs,
                       pheno = all_df_shared_factors_mat_rankZ,
                       kinship = shared.kinship,
                       addcovar = shared.covar)

# get peaks above lod 5
all_df_shared_peaks <- find_peaks( all_df_scans, map = gmap, threshold = 5)
# interpolate physical location of peaks 
all_df_shared_peaks <- all_df_shared_peaks %>%
  mutate(phenotype=lodcolumn) %>%
  mutate( peak_chr = chr,
          peak_cM = pos,
          factor_num = as.numeric(str_sub(lodcolumn, 7)) )%>%
  interp_bp(.) #add bp location for peaks

# add markers before & after the QTL peaks for getting effects later
# Get the bounding markers for each QTL peak
# i.e. markers on the 69k grid that are up- and downstream of the peak
query <- all_df_shared_peaks %>% 
  dplyr::select(peak_chr, interp_bp_peak) %>%
  dplyr::rename(chrom=peak_chr, start=interp_bp_peak) %>% mutate(end=start) %>%
  GenomicRanges::GRanges()
subject <- select(map_dat2, chrom, pos_bp) %>% 
  dplyr::rename(start=pos_bp) %>%
  mutate(end=start) %>% 
  GenomicRanges::GRanges()   # length 69,005

all_df_shared_peaks$before <- map_dat2$marker[follow(query, subject)]
all_df_shared_peaks$after <- map_dat2$marker[precede(query, subject)]

# loading the permutation results for significance thresholds
load(here("../pQTL_website/_data/MOFA_all_df_thres_rankZ.RData"))


# get the list of significant MOFA QTL
summary(all_df_shared_thres) %>%
  as_tibble( rownames = "Significance level") %>%
  pivot_longer(names_to = "Factor", values_to="thres", cols = 2:24) %>%
  left_join( all_df_shared_peaks, by = c("Factor"= "phenotype")) %>%
  mutate( phenotype = Factor) %>%
  mutate( lod = round(lod,1), thres = round(thres, 1))  %>%
  filter( lod >= thres  ) %>%
  filter( !(chr =="X" & lod <= thres+5) ) -> all_df_shared_peaks_sig

all_df_shared_peaks_sig %>% 
  mutate( peak_pos_Mbp = round((interp_bp_peak/1e06),2),
          lod = round(lod,1),
          factor_num = paste("Factor", factor_num),
          qtl_loc = paste0(chr, ": ",peak_pos_Mbp)) %>% 
  select( `MOFA Factor` = factor_num,
          `QTL Peak location (Chr: Mbp)` = qtl_loc,
          `LOD score` = lod) %>% 
  mutate(
    Type = case_when( 
      `MOFA Factor` %in% c("Factor 3","Factor 13")~"caQTL, eQTL",
      `MOFA Factor` == "Factor 4"~"eQTL",
      `MOFA Factor` %in% c("Factor 12","Factor 16")~"New",
      `MOFA Factor` == "Factor 14" & `QTL Peak location (Chr: Mbp)`=="16: 65.58"~"New",
      `MOFA Factor` == "Factor 14" & `QTL Peak location (Chr: Mbp)`=="4: 147.89"~"caQTL, eQTL",
      )
  ) %>% 
  create_dt()

Data shown in Figure 5B.

Figure 5C: MOFA Factor 3 QTL on chromosome 15

# get allele effects for MOFA Factor 3 on chr 15
factor3_chr15_qtl <- all_df_shared_peaks_sig %>%
  filter( factor_num ==3, chr == 15)

# run effects scan using qtl2::scan1blup & get the mean of allele effects between QTL borders (markers before & after)
factor3_chr15_qtl_effs_scan <- scan1blup(
  genoprobs =shared.probs[,15],
  pheno = all_df_shared_factors_mat_rankZ[,"Factor3", drop = FALSE],
  kinship = shared.kinship[[15]],
  addcovar = shared.covar)
factor3_chr15_qtl_effs <- colMeans(factor3_chr15_qtl_effs_scan[c(factor3_chr15_qtl$before, factor3_chr15_qtl$after), LETTERS[1:8]]) 

# # get lod scores + effects for all the proteins at the same locus
pQTL_at_chr15_peak <- c()
mofa_factor <- "Factor3"
mofa_peak_chr <- factor3_chr15_qtl$chr
mofa_markers <- c(factor3_chr15_qtl$before, factor3_chr15_qtl$after)
probs2_markers <- subset_probs( probs.esc_prot,
                                this_chrom =mofa_peak_chr,
                                this_markers = mofa_markers )

for( i in 1:nrow(all.prots)){
    probs2_markers <- subset_probs( probs.esc_prot,
                                    this_chrom = mofa_peak_chr,
                                    this_markers = mofa_markers )
    pQTL_eff <- scan1blup(genoprobs = probs2_markers,
                          pheno = exprZ.esc_prot[, all.prots$protein_id[i],drop=FALSE],
                          kinship_loco.esc_prot[[mofa_peak_chr]],
                          addcovar = covar.esc_prot
    )
    pQTL_mean_eff <- colMeans(pQTL_eff[,LETTERS[1:8]])
    names(pQTL_mean_eff) <- paste0(LETTERS[1:8],".esc_prot")

    pQTL_at_chr15_peak[[i]] <- c( protein_id = all.prots$protein_id[i],
                                  peak_chr = mofa_peak_chr,
                                  Factor = mofa_factor,
                                  pQTL_mean_eff)

  }

# pQTL
# convert to matrix
pQTL_at_chr15_peak_effs <- pQTL_at_chr15_peak %>% 
  do.call(rbind,.) %>% 
  as_tibble() %>% 
  mutate_at(., c(paste0(LETTERS[1:8],".esc_prot")), as.numeric )
pQTL_at_chr15_peak_eff_mat <- pQTL_at_chr15_peak_effs%>% 
  select(c(protein_id,paste0(LETTERS[1:8],".esc_prot")) ) %>% 
  column_to_rownames("protein_id") %>%
  as.matrix() %>% 
  t()
rownames(pQTL_at_chr15_peak_eff_mat) <- LETTERS[1:8]
factor3_chr15_qtl_effs_mat <- factor3_chr15_qtl_effs %>% 
  as.matrix()
# get correlations
factor3_pqtl_effs_corr <- cor(pQTL_at_chr15_peak_eff_mat, factor3_chr15_qtl_effs_mat)
factor3_pqtl_effs_wcorr <- factor3_pqtl_effs_corr %>% 
  as_tibble( rownames = "protein_id") %>% 
  rename( "pqtl_cor"= V1) 

# add pQTL lod for each protein
# get pQTL lods from scans + add ids
# load scans
load("/projects/munger-lab/projects/DO_mESC/proteomics/pqtl_mapping_SA/DO195_mESC_pQTL_scans_noPoly_v2.RData") # esc.prot.scans
factor3_pqtl_lods <- apply( esc.prot.scans[c(factor3_chr15_qtl$before, factor3_chr15_qtl$after),],2,max ) %>% 
  as_tibble( rownames = "protein_id") %>% 
  rename( pQTL_lod = value)
rm(esc.prot.scans)

# merge all into a data frame for plotting
Figure5C_data <- factor3_pqtl_effs_wcorr %>% 
  left_join( factor3_pqtl_lods) %>% 
  left_join( all_df_shared_weights %>% 
               filter( factor == "Factor3", feature %in% all.prots$protein_id) %>% 
               select( protein_id = feature, protein_weight = value)
  ) %>%
  # add protein details
  left_join( ., all.prots2 %>% 
               select( protein_id, gene_chr, midpoint, mgi_symbol)) %>% 
  # add mofa qtl details
  cbind( ., factor3_chr15_qtl %>% 
           select( factor3_qtl_pos =interp_bp_peak , factor3_qtl_peak_chr = peak_chr)
         ) %>% 
  filter(!is.na(pqtl_cor)) %>% 
  filter(  !(factor3_qtl_peak_chr == gene_chr &
              abs( midpoint - as.numeric(factor3_qtl_pos)) < 10e06 ) ) %>% 
  select(
    `Protein ID` = protein_id,
    `MGI symbol` = mgi_symbol,
    `Correlation between QTL effects` = pqtl_cor,
    `LOD score` = pQTL_lod,
    `Protein weight` = protein_weight
  )

Figure5C_data %>% 
  filter( `LOD score` < 7.5, 
          abs(`Protein weight`) >0.1 ) %>% 
  ggplot()+
  aes(
    x = `Protein weight`,
    y =  `LOD score`,
    col = `Correlation between QTL effects`,
    text = paste("Gene: ", `MGI symbol`) 
  )+
  geom_point(size = 3, alpha = 0.6)+
  theme_pubclean(base_size = 18)+
  scale_color_gradient2( limits = c(-1,1))+
  xlab("Protein weight")+
  ylab("LOD score")+
  labs( col = "Correlation between\nQTL effects") +
  xlim(-0.75,0.4)+
  ylim(0,15)+
  theme(legend.text = element_text(angle = 30), 
        legend.position = "top",
        legend.title = element_text(vjust = 1.3))+
  # add the signifincatn qtl in gray
  geom_point( data = filter(Figure5C_data, `LOD score` > 7.5, abs(`Protein weight`)>0.1),
              col = "gray", size =3, alpha = 0.6)-> Figure5C

#ggplotly(Figure5C,  width = 700, height = 400)
 
Figure5C 

Figure 5C: For all expressed proteins, the pQTL LOD score calculated at the Chr 15 QTL peak is plotted on the y-axis relative to the protein’s contribution (factor weight) to MOFA Factor 3 on the x-axis. Proteins with absolute factor weights less than 0.1 were filtered. Correlation between allele effects at the Chr 15 pQTL for individual proteins to allele effects of the Factor 3 QTL. Individual genes that mapped with a significant QTL (LOD > 7.5) are colored gray, and highlight that many proteins contribute substantially to Factor 3 and show high agreement in allele effects at the Chr 15 QTL (dark red and blue), despite not mapping individually with a significant QTL at that locus.

Figure5C_data %>% 
  mutate_if( is.numeric, round, 2) %>% 
  create_dt()

Data plotted in Figure 5C.

Figure 5D: MOFA Factor 4 QTL on chromosome 10

# get allele effects for MOFA Factor 4 on chr 10
factor4_chr10_qtl <- all_df_shared_peaks_sig %>%
  filter( factor_num == 4, chr == 10)

# run effects scan using qtl2::scan1blup & get the mean of allele effects between QTL borders (markers before & after)
factor4_chr10_qtl_effs_scan <- scan1blup(
  genoprobs =shared.probs[,10],
  pheno = all_df_shared_factors_mat_rankZ[,"Factor4", drop = FALSE],
  kinship = shared.kinship[[10]],
  addcovar = shared.covar)
factor4_chr10_qtl_effs <- colMeans(factor4_chr10_qtl_effs_scan[c(factor4_chr10_qtl$before, factor4_chr10_qtl$after), LETTERS[1:8]])

# get lod scores + effects for all the proteins at the same locus
eQTL_at_chr10_peak <- c()
mofa_factor <- "Factor4"
mofa_peak_chr <- factor4_chr10_qtl$chr  
mofa_markers <- c(factor4_chr10_qtl$before, factor4_chr10_qtl$after)
probs2_markers <- subset_probs( probs.esc_rna, 
                                this_chrom =mofa_peak_chr, 
                                this_markers = mofa_markers )

for( i in 1:nrow(all.genes)){
    probs2_markers <- subset_probs( probs.esc_rna, 
                                    this_chrom = mofa_peak_chr, 
                                    this_markers = mofa_markers )
  
    eQTL_eff <- scan1blup(genoprobs = probs2_markers,
                          pheno = exprZ.esc_rna[, all.genes$ensembl_gene_id[i],drop=FALSE],
                          kinship_loco.esc_rna[[mofa_peak_chr]],
                          addcovar = covar.esc_rna
    )
    eQTL_mean_eff <- colMeans(eQTL_eff[,LETTERS[1:8]])
    names(eQTL_mean_eff) <- paste0(LETTERS[1:8],".esc_rna")
    
    eQTL_at_chr10_peak[[i]] <- c( ensembl_gene_id = all.genes$ensembl_gene_id[i],
                                  peak_chr = mofa_peak_chr,
                                  Factor = mofa_factor,
                                  eQTL_mean_eff)
    
  }

# eQTL
# convert to matrix
eQTL_at_chr10_peak_effs <- eQTL_at_chr10_peak %>% 
  do.call(rbind,.) %>% 
  as_tibble() %>% 
  mutate_at(., c(paste0(LETTERS[1:8],".esc_rna")), as.numeric )
eQTL_at_chr10_peak_effs_mat <- eQTL_at_chr10_peak_effs%>% 
  select(c(ensembl_gene_id,paste0(LETTERS[1:8],".esc_rna")) ) %>% 
  column_to_rownames("ensembl_gene_id") %>%
  as.matrix() %>% 
  t()
rownames(eQTL_at_chr10_peak_effs_mat) <- LETTERS[1:8]
factor4_chr10_qtl_effs_mat <- factor4_chr10_qtl_effs %>% 
  as.matrix()
# get correlations
factor4_eqtl_effs_corr <- cor(eQTL_at_chr10_peak_effs_mat, factor4_chr10_qtl_effs_mat)
factor4_eqtl_effs_wcorr <- factor4_eqtl_effs_corr %>% 
  as_tibble( rownames = "ensembl_gene_id") %>% 
  rename( "eqtl_cor"= V1) 

# add eQTL lod
# get eQTL lods from scans + add ids
# load scans
load("/projects/munger-lab/projects/DO_mESC/rna_seq/qtl_mapping/total_gene_expression/eqtl_grid69k_pe/DO185_mESC_paired_eQTL_scans.RData") # esc.rna.scans
factor4_eqtl_lods <- apply( esc.rna.scans[c(factor4_chr10_qtl$before, factor4_chr10_qtl$after),],2,max ) %>% 
  as_tibble( rownames = "ensembl_gene_id") %>% 
  rename( eQTL_lod = value)
rm(esc.rna.scans)

# merge all into a data frame
Figure5D_data <- factor4_eqtl_effs_wcorr %>% 
  left_join( factor4_eqtl_lods) %>% 
  left_join( all_df_shared_weights %>% 
               filter( factor == "Factor4", feature %in% all.genes$ensembl_gene_id) %>% 
               select( ensembl_gene_id = feature, transcript_weight = value)
  ) %>% 
  # add gene details
  left_join( ., all.genes2 %>% 
               select( ensembl_gene_id, midpoint, gene_chr, mgi_symbol)) %>% 
  # add mofa qtl details
  cbind(., factor4_chr10_qtl %>% 
          select( factor4_qtl_pos =interp_bp_peak , factor4_qtl_peak_chr = peak_chr)
  ) %>% 
  filter(!is.na(eqtl_cor)) %>% 
  filter(  !(factor4_qtl_peak_chr == gene_chr &
              abs( midpoint - as.numeric(factor4_qtl_pos)) < 10e06 ) ) %>% 
  select(
    `Gene ID` = ensembl_gene_id,
    `MGI symbol` = mgi_symbol,
    `Correlation between QTL effects` = eqtl_cor,
    `LOD score` = eQTL_lod,
    `Transcript weight` = transcript_weight
  )

Figure5D_data %>% 
  filter( `LOD score` < 7.5, 
          abs(`Transcript weight`) >0.1 ) %>% 
  ggplot()+
  aes(
    x = `Transcript weight`,
    y =  `LOD score`,
    col = `Correlation between QTL effects`,
    text = paste("Gene: ", `MGI symbol`) 
  )+
  geom_point(size = 3, alpha = 0.6)+
  theme_pubclean(base_size = 18)+
  scale_color_gradient2( limits = c(-1,1))+
  xlab("Transcript weight")+
  ylab("LOD score")+
  labs( col = "Correlation between\nQTL effects") +
  ylim(0,15)+
  xlim(-0.5, 1)+
  theme(legend.text = element_text(angle = 30), 
        legend.position = "top",
        legend.title = element_text(vjust = 1.3))+
  # add the significant qtl in gray
  geom_point( data = filter(Figure5D_data, `LOD score` > 7.5, abs(`Transcript weight`)>0.1),
              col = "gray", size =3, alpha = 0.6) -> Figure5D

#ggplotly(Figure5D, width = 700, height = 400)
Figure5D

Figure 5D: For all expressed transcripts, the eQTL LOD score at the Chr 10 QTL peak is plotted on the y-axis relative to that transcript’s contribution to Factor 4 on the x-axis. Again, transcripts with absolute factor weights less than 0.1 were filtered, and individual points are as described in panel C. Many transcripts contribute to Factor 4 and have correlated allele effects at the Chr 10 QTL, despite individually failing to map with a significant Chr 10 eQTL.

Figure5D_data %>% 
  mutate_if( is.numeric, round, 2) %>% 
  create_dt()

Data plotted in Figure 5D.

Figure 5E: Mediation of Factor 4 QTL on chromosome 10

# prep factor 4 qtl for mediation
factor4_chr10_qtl <- factor4_chr10_qtl %>% 
  mutate(pos_cM = as.double(pos)) %>%
  left_join(select(map_dat2, -chr), by = c("pos_cM")) %>%
  filter(!is.na(marker)) %>%
  mutate(phenotype = lodcolumn)

# factor 4 qtl scan
scan1_factor4 <- scan1(genoprobs = shared.probs,
                 pheno =  all_df_shared_factors_mat_rankZ[,"Factor4",drop = FALSE],
                 kinship = shared.kinship,
                 addcovar = shared.covar)

# factor 4 mediation
# phenotype to be mediated  
target <- all_df_shared_factors_mat_rankZ[, "Factor4", drop = FALSE]
# prep ESC transcript mediators & annotations
meds.rna <- tibble(ensembl_gene_id = colnames(expr.esc_rna)) %>%
  left_join(all.genes) %>%
  mutate( id = ensembl_gene_id) %>%
  mutate(chrom = ifelse(gene_chr == "MT", "M", gene_chr)) %>%
  filter(!is.na(chrom))
mediator.rna <- exprZ.esc_rna[threeway.shared.samples$sampleid,
  meds.rna$ensembl_gene_id,
  drop = FALSE
] %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  left_join(., select(threeway.shared.samples, sampleid, top_muga), by = c("rowname" = "sampleid")) %>%
  column_to_rownames("top_muga") %>%
  select(-rowname) %>%
  as.matrix()
# prep ESC protein mediators & annotations
meds.prot <- tibble(protein_id = colnames(exprZ.esc_prot[,all.prots$protein_id])) %>%
  left_join(all.prots) %>%
  mutate( id = protein_id) %>%
  mutate(chrom = ifelse(gene_chr == "MT", "M", gene_chr)) %>%
  mutate(end = gene_end, start = gene_start, id = protein_id)
mediator.prot <- exprZ.esc_prot[threeway.shared.samples$sampleid,
  meds.prot$protein_id,
  drop = FALSE
] %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  left_join(., select(threeway.shared.samples, sampleid, top_muga), by = c("rowname" = "sampleid")) %>%
  column_to_rownames("top_muga") %>%
  select(-rowname) %>%
  as.matrix()
# run mediation using a wrapper function
get_meds <- function(target, meds, mediator, peaks, probs,covar, z_thres = -4,  pos_thres = 10 ){
  med.scan.all <- c()
  samples <- rownames(mediator)
  for( i in 1:nrow(peaks)){
    #print( paste0("i is ", i))
    
    marker    <- map_dat2 %>% filter(pos_cM == peaks$pos_cM[i])
    annot     <- meds %>% mutate(chr=gene_chr,pos=abs(gene_end+gene_start)/2)
    geno      <- pull_genoprobpos(probs,marker$marker)
    geno      <- geno[samples,]
    
    if( !is.null(covar)){ 
      covar <- covar[samples,,drop=FALSE] 
    }
    
    med.scan <- mediation.scan(target= target[samples,peaks$phenotype[i], drop=FALSE],
                               mediator = mediator,
                               annotation = annot,
                               covar =  covar,
                               qtl.geno = geno,
                               verbose = FALSE,
                               method     = "double-lod-diff") 
    
    med.scan <- med.scan %>% 
      mutate(phenotype= peaks$phenotype[i], 
             peak_chr = peaks$peak_chr[i], 
             peak_lod = peaks$lod[i],
             peak_pos_Mbp = peaks$pos_bp/1e06,
             med_chr  = chr) %>% 
      mutate(scaled_LOD = scale(LOD), 
             middle = (gene_end+gene_start)/2e06) %>%
      filter( (#scaled_LOD < z_thres & 
                 peak_chr ==  med_chr & 
                 abs(middle - peak_pos_Mbp) <= pos_thres) &
                 (LOD < peak_lod))
    
    med.scan.all <- rbind(med.scan.all,med.scan)
  }
  return(med.scan.all)
}
factor4_rna_meds <-  get_meds(
  target = target,
  meds = meds.rna,
  mediator = mediator.rna,
  peaks = factor4_chr10_qtl,
  probs = shared.probs,
  covar = shared.covar
)
factor4_prot_meds <-  get_meds(
  target = target,
  meds = meds.prot,
  mediator = mediator.prot,
  peaks = factor4_chr10_qtl,
  probs = shared.probs,
  covar = shared.covar
)

factor4_meds <- factor4_rna_meds %>% 
  mutate( type = "RNA") %>% 
  rbind( factor4_prot_meds %>% 
           mutate( type = "protein")) %>% 
  mutate(
    `Phenotype` = phenotype, 
    `Mediator MGI symbol` = mgi_symbol,
    `QTL chromosome` = peak_chr, 
    `Mediator chromosome` = chr,
    `QTL LOD` = peak_lod,
    `Mediated LOD` = LOD,
    `QTL position (Mbp)`= peak_pos_Mbp, 
    `Mediator midpoint (Mbp)` = middle
  )

scan1_factor4_df <- scan1_factor4 %>% 
  as.data.frame( ) %>% 
  rename( `Factor 4 QTL` = "Factor4") %>% 
  mutate( marker = dimnames(scan1_factor4)[[1]]) %>% 
  left_join(map_dat2) 
scan1_factor4_df$cumsum_pos_bp <- scan1_factor4_df$pos_bp + chrom_lens_offset[scan1_factor4_df$chr]
scan1_factor4_df <- scan1_factor4_df %>% 
  mutate( `Offestted chromosome coordinates` = cumsum_pos_bp,
          marker,
          Chromosome = chr,
          `Coordinates (cM)` = pos_cM, 
          `Coordinates (bp)` = pos_bp) 

Figure5E_data_qtl_scan <- scan1_factor4_df %>% 
  select( `Factor 4 QTL`, 
          `Offestted chromosome coordinates` = cumsum_pos_bp,
          marker,
          Chromosome = chr,
          `Coordinates (cM)` = pos_cM, 
          `Coordinates (bp)` = pos_bp)

Figure5E_data_mediation <- factor4_meds %>% 
  select(
    `Phenotype` = phenotype, 
    `Mediator MGI symbol` = mgi_symbol,
    `QTL chromosome` = peak_chr, 
    `Mediator chromosome` = chr,
    `QTL LOD` = peak_lod,
    `Mediated LOD` = LOD,
    `QTL position (Mbp)`= peak_pos_Mbp, 
    `Mediator midpoint (Mbp)` = middle,
    type
  )

# graphical prep for chromosome locations
chroms <- c(as.character(1:19), "X")
chrom_lens <- c( 195431559, 182107670, 160017104, 156496071, 151833620, 149721874, 145434693, 129399468, 124582650, 130685419, 122078650, 120120622 ,120387272, 124867725, 104015452, 98180002, 94984432, 90672596, 61417310 , 171028300)
names(chrom_lens) <- chroms
chrom_lens_offset <- cumsum(chrom_lens) - chrom_lens
chrom_lens_midpt <- chrom_lens_offset + chrom_lens / 2

# lod plot 
Figure5E_data_qtl_scan %>% 
  ggplot()+
    aes( 
      x= `Offestted chromosome coordinates`,
      y = `Factor 4 QTL`
      )+
  geom_line( size = 1, col = "dark gray",alpha = 0.8)+
  theme_pubclean( base_size = 18)+
  xlab("Chromosome")+
  ylab( "LOD score")+
  scale_x_discrete( name = "Chromosome",
                    limits = chrom_lens_midpt, 
                    labels = names(chrom_lens), 
                    expand = expansion( mult = 0.05)) -> factor4_lod_plot

# overlay mediation
# highlights
factor4_meds_min <- Figure5E_data_mediation %>% 
  slice_min( `Mediated LOD`)
factor4_meds_duxf <- Figure5E_data_mediation %>% 
  filter(`Mediator MGI symbol` == "Duxf3")

factor4_lod_plot+
  geom_point( data = Figure5E_data_mediation, 
              aes( x = `Mediator midpoint (Mbp)`*1e06+chrom_lens_offset[10],
                   y = `Mediated LOD`),
              col = qtl.colors[["rna"]],
              size = 4,
              alpha = 0.9)+
  annotate( "label", 
            y = factor4_meds_min$`Mediated LOD`+0.1, 
            x = factor4_meds_min$`Mediator midpoint (Mbp)`*1e06+chrom_lens_offset[10]-1.9e08 ,
            label = factor4_meds_min$`Mediator MGI symbol`, 
            size = 6, 
            col = qtl.colors[["rna"]],
            fontface = "italic")+
  annotate("label",
           y = factor4_meds_duxf$`Mediated LOD`+0.1, 
            x =factor4_meds_duxf$`Mediator midpoint (Mbp)`*1e06+chrom_lens_offset[10]-1.3e08 ,
            label = factor4_meds_duxf$`Mediator MGI symbol`, 
            size = 6, 
            col = qtl.colors[["rna"]],
            fontface = "italic")+
  ggtitle("MOFA Factor 4 QTL")-> factor4_med_plot

factor4_med_plot

Figure 5E:Genome-wide LOD scores obtained from the Factor 4 QTL scan is plotted with mediation results overlaid. Mediation with Gm20625 transcript abundance causes the largest decrease in QTL LOD score. Duxf3 expression was previously identified as a strong candidate mediator for the eQTL hotspot in this region (Skelly et al., 2020), but performs poorly as a mediator of the Factor 4 QTL compared to Gm20625. Both genes are highlighted in green next to their corresponding LOD score drop.

QTL scan used in plotting Figure 5E can be downloaded below.

list(Figure5E_data_qtl_scan) %>% 
  downloadthis::download_this(    output_name = "Figure5E data",
    output_extension = ".xlsx",
    button_label = "Download Figure 5E data as xlsx",
    button_type = "primary",
    has_icon = TRUE,
    icon = "fa fa-save"
    )

Download Figure 5E data as xlsx

Figure5E_data_mediation %>% 
  mutate_if(is.numeric, round, 2) %>% 
  arrange(`Mediated LOD`) %>% 
  create_dt()

Mediation results plotted in Figure 5E.

---

A work by Selcan Aydin

selcan.aydin@jax.org